tidy_custom.polr <- function(x, ...) {
  s <- coeftest(x)
  out <- data.frame(
    term = row.names(s),
    p.value = s[, "Pr(>|t|)"])
  out
}

LDV_table_function <- function(model,
                               data_frames,
                               outcomes,
                               treatments,
                               trust_covariates,
                               interaction_effects = NULL,
                               interaction_effects_name = NULL,
                               fixed_effects_name,
                               table_name,
                               return_object = FALSE) {
  
  n_models <- length(data_frames)
  
  if (length(interaction_effects) == 0) {
    
    models <- lapply(1:n_models, function(i)
      
      if (model == "ologit") {
        if (identical(data_frames[[i]], diplomacy_long)) {
          polr(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], " + factor(", trust_covariates[i], ")", "*factor(country_exp)", " + factor(", fixed_effects_name, ")")),
            Hess = TRUE,
            data = data_frames[[i]])
        } else {
          polr(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], " + factor(", trust_covariates[i], ")" , " + factor(", fixed_effects_name, ")")
          ),
          Hess = TRUE,
          data = data_frames[[i]]) 
        }
      } else if (model == "logit") {
          if (identical(data_frames[[i]], diplomacy_long)) {
            glm(as.formula(
              paste0("factor(", outcomes[i], ") ~ ", treatments[i], " + factor(", trust_covariates[i], ")", "*factor(country_exp)", " + factor(", fixed_effects_name, ")")),
              family = "binomial",
              data = data_frames[[i]])
          } else {
            glm(as.formula(
              paste0("factor(", outcomes[i], ") ~ ", treatments[i], " + factor(", trust_covariates[i], ")" , " + factor(", fixed_effects_name, ")")
            ),
            family = "binomial",
            data = data_frames[[i]]) 
          }
      }

  )
    
  } else if (length(interaction_effects) == 1) {
    
    models <- lapply(1:n_models, function(i)
      
      if (model == "ologit") {
        if (identical(data_frames[[i]], diplomacy_long)) {
          polr(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], "*", interaction_effects[[1]][i],
                     " + factor(", trust_covariates[i], ")", "*factor(country_exp)", " + factor(", fixed_effects_name, ")")
          ),
          Hess = TRUE,
          data = data_frames[[i]])
        } else {
          polr(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], "*", interaction_effects[[1]][i], " + factor(", trust_covariates[i], ") + ", "factor(", fixed_effects_name, ")")
          ),
          Hess = TRUE,
          data = data_frames[[i]])
        }
      } else if (model == "logit") {
        if (identical(data_frames[[i]], diplomacy_long)) {
          glm(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], "*", interaction_effects[[1]][i],
                   " + factor(", trust_covariates[i], ")", "*factor(country_exp)", " + factor(", fixed_effects_name, ")")
          ),
          family = "binomial",
          data = data_frames[[i]])
        } else {
          glm(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], "*", interaction_effects[[1]][i], " + factor(", trust_covariates[i], ") + ", "factor(", fixed_effects_name, ")")
          ),
          family = "binomial",
          data = data_frames[[i]])
        }
      }
    )
    
  } else if (length(interaction_effects) == 2) {
    
    models <- lapply(1:n_models, function(i)
      
      if (model == "ologit") {
        if (identical(data_frames[[i]], diplomacy_long)) {
          polr(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i], " + factor(", trust_covariates[i], ")", "*factor(country_exp)", " + factor(", fixed_effects_name, ")")
          ),
          Hess = TRUE,
          data = data_frames[[i]])
        } else {
          polr(as.formula(
            paste0("factor(", outcomes[i], ") ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i], " + factor(", trust_covariates[i], ")", " + factor(", fixed_effects_name, ")")
          ),
          Hess = TRUE,
          data = data_frames[[i]])
        }
      } else if (model == "logit") {
        if (identical(data_frames[[i]], diplomacy_long)) {
          glm(as.formula(
            paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i], " + factor(", trust_covariates[i], ")", "*factor(country_exp)", " + factor(", fixed_effects_name, ")")
          ),
          family = "binomial",
          data = data_frames[[i]])
        } else {
          glm(as.formula(
            paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i], " + factor(", trust_covariates[i], ")", " + factor(", fixed_effects_name, ")")
          ),
          family = "binomial",
          data = data_frames[[i]])
        }
      }
    ) 
  }
  
  names(models) <- names(outcomes)
  
  models_frame <- lapply(1:n_models, function(i)
    model.frame(models[[i]])
  )
  
  for (i in 1:length(models)) {
    names(models_frame[[i]]) <- str_replace_all(names(models_frame[[i]]), "factor\\(", "")
    names(models_frame[[i]]) <- str_replace_all(names(models_frame[[i]]), "\\)", "")
  }
  
  outcome_stats <- lapply(1:n_models, function(i)
    filter(models_frame[[i]], get(treatments[i]) == 0)
  )
  
  outcome_stats <- lapply(1:n_models, function(i)
    outcome_stats[[i]] %>%
      mutate(across(unname(outcomes[i]), ~ as.numeric(as.character(.)))) %>%
      dplyr::summarize(Mean = round(mean(get(outcomes[i]), na.rm = TRUE),2),
               SD = round(sd(get(outcomes[i]), na.rm = TRUE),2)))
  
  outcome_stats <- lapply(1:n_models, function(i)
    outcome_stats[[i]] %>%
      mutate(Range = paste(sort(unique(models_frame[[i]][,outcomes[i]])), collapse = ",")) %>%
      mutate(Range = paste0("\\{", Range, "\\}"))
  )
  
  names(outcome_stats) <- names(outcomes)
  
  treat_mean <- as.numeric(sapply(1:n_models, function(i)
    models_frame[[i]] %>%
      dplyr::summarize(mean = round(mean(get(treatments[i]), na.rm = TRUE),2))
  ))
  
  if (length(interaction_effects) != 0) {
    
    interaction_stats <-
      lapply(1:length(interaction_effects), function (i)
        lapply(1:n_models, function(j)
          models_frame[[j]] %>%
            dplyr::summarize(
              Mean = round(mean(get(
                interaction_effects[[i]][j]
              ), na.rm = TRUE), 2),
              SD = round(sd(get(
                interaction_effects[[i]][j]
              ), na.rm = TRUE), 2),
              Min = round(min(get(
                interaction_effects[[i]][j]
              ), na.rm = TRUE), 2),
              Max = round(max(get(
                interaction_effects[[i]][j]
              ), na.rm = TRUE), 2)
            )))
    
  }
  
  if (length(interaction_effects) == 0) {
    
    extra_rows <- rbind(
      sapply(outcome_stats, function(x) x[,'Range']),
      sapply(outcome_stats, function(x) x[,'Mean']),
      sapply(outcome_stats, function(x) x[,'SD']),
      treat_mean
    )
    
    extra_rows <- as_tibble(extra_rows) %>%
      mutate(term = c("Outcome range",
                      "Control outcome mean",
                      "Control outcome std. dev.",
                      paste0(unique(names(treatments)), " mean"))) %>%
      relocate(term)
    
    attr(extra_rows, 'position') <- 3:6
    
    treatment_names <- treatments[!duplicated(treatments)]
    coef_map_list <- setNames(names(treatment_names), treatment_names)
    
  } else if (length(interaction_effects) == 1) {
    
    extra_rows <- rbind(
      paste0("[", sapply(interaction_stats[[1]], function(x) x[,"Min"]),
             ",", sapply(interaction_stats[[1]], function(x) x[,"Max"]), "]"),
      sapply(interaction_stats[[1]], function(x) x[,'Mean']),
      sapply(interaction_stats[[1]], function(x) x[,'SD']),
      sapply(outcome_stats, function(x) x[,'Range']),
      sapply(outcome_stats, function(x) x[,'Mean']),
      sapply(outcome_stats, function(x) x[,'SD']),
      treat_mean
    )
    
    extra_rows <- as_tibble(extra_rows) %>%
      mutate(term = c(paste0(interaction_effects_name[[1]], " range"),
                      paste0(interaction_effects_name[[1]], " mean"),
                      paste0(interaction_effects_name[[1]], " std. dev."),
                      "Outcome range",
                      "Control outcome mean",
                      "Control outcome std. dev.",
                      paste0(unique(names(treatments)), " mean"))) %>%
      relocate(term)
    
    coef_map_list <- c(names(treatments[!duplicated(treatments)]),
                       rep(paste0(names(treatments[!duplicated(treatments)]), " $\\times$ ", interaction_effects_name[[1]]), n_models))
    
    names(coef_map_list) <- c(treatments[!duplicated(treatments)],
                              paste0(treatments[!duplicated(treatments)], ":", interaction_effects[[1]]))
    
    coef_map_list <- coef_map_list[!duplicated(names(coef_map_list))]
    
    attr(extra_rows, 'position') <- 5:12
    
  } else if (length(interaction_effects) == 2) {
    
    extra_rows <- rbind(
      paste0("[", sapply(interaction_stats[[1]], function(x) x[,"Min"]),
             ",", sapply(interaction_stats[[1]], function(x) x[,"Max"]), "]"),
      sapply(interaction_stats[[1]], function(x) x[,'Mean']),
      sapply(interaction_stats[[1]], function(x) x[,'SD']),
      paste0("[", sapply(interaction_stats[[2]], function(x) x[,"Min"]),
             ",", sapply(interaction_stats[[2]], function(x) x[,"Max"]), "]"),
      sapply(interaction_stats[[2]], function(x) x[,'Mean']),
      sapply(interaction_stats[[2]], function(x) x[,'SD']),
      sapply(outcome_stats, function(x) x[,'Range']),
      sapply(outcome_stats, function(x) x[,'Mean']),
      sapply(outcome_stats, function(x) x[,'SD']),
      treat_mean
    )
    
    extra_rows <- as_tibble(extra_rows) %>%
      mutate(term = c(paste0(interaction_effects_name[[1]], " range"),
                      paste0(interaction_effects_name[[1]], " mean"),
                      paste0(interaction_effects_name[[1]], " std. dev."),
                      paste0(interaction_effects_name[[2]], " range"),
                      paste0(interaction_effects_name[[2]], " mean"),
                      paste0(interaction_effects_name[[2]], " std. dev."),
                      "Outcome range",
                      "Control outcome mean",
                      "Control outcome std. dev.",
                      paste0(unique(names(treatments)), " mean"))) %>%
      relocate(term)
    
    coef_map_list <- c(names(treatments[!duplicated(treatments)]),
                       rep(paste0(names(treatments[!duplicated(treatments)]), " $\\times$ ", interaction_effects_name[[1]]), n_models),
                       rep(paste0(names(treatments[!duplicated(treatments)]), " $\\times$ ", interaction_effects_name[[2]]), n_models))
    
    names(coef_map_list) <- c(treatments[!duplicated(treatments)],
                              paste0(treatments[!duplicated(treatments)], ":", interaction_effects[[1]]),
                              paste0(treatments[!duplicated(treatments)], ":", interaction_effects[[2]]))
    
    coef_map_list <- coef_map_list[!duplicated(names(coef_map_list))]
    
    attr(extra_rows, 'position') <- 7:16
    
  }
  
  extra_rows <- rbind(extra_rows,
                      c("Number of blocks", rep("Number of blocks not provided for LDV models", n_models)))
  
  attr(extra_rows, 'position') <- c(attr(extra_rows, 'position'), max(attr(extra_rows, 'position'))+1)
  
  f <- function(x) format(round(x, 2), big.mark=",")
  gof_mapping <- list(
    list("raw" = "r.squared", "clean" = "R$^2$", "fmt" = f),
    list("raw" = "nobs", "clean" = "Observations", "fmt" = f)
  )
  
  if (model == "logit") {
    
  models <- lapply(1:n_models, function(i)
    if (identical(data_frames[[i]], diplomacy_long)) {
      coeftest(models[[i]], vcovCL(models[[i]], cluster = data_frames[[i]]$mergeid))
    } else {
      coeftest(models[[i]], vcovHC(models[[i]], type = "HC1"))
    }
  )
  
  names(models) <- names(outcomes)
    
  }
  
  modelsummary(models,
               coef_map = coef_map_list,
               add_rows = extra_rows,
               gof_omit = "se_type|Adj.",
               gof_map = gof_mapping,
               stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
               escape = FALSE,
               output = paste0("Tables/", table_name, ".tex"))
  
  if (return_object == TRUE) {
    
    table_object <- list(models = models,
                         extra_rows = extra_rows)
    
    return(table_object)
    
  } 
  
}